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Abstract. This paper considers a generalization of Gaussian random field 
with covariance function of Whittle-Matern family. Such a random field can 
be obtained as the solution to the fractional stochastic differential equation 
with two fractional orders. Asymptotic properties of the covariance functions 
belonging to this generalized Whittle-Matern family are studied, which are 
used to deduce the sample path properties of the random field. The Whittle- 
Matern field has been widely used in modeling geostatistical data such as sea 
beam data, wind speed, field temperature and soil data. In this article we 
show that generalized Whittle-Matern field provides a more flexible model for 
wind speed data. 



1. Introduction 

Random fields play an important role in geostatistics, which deals with problems 
stretching from the resource evaluation such as the estimation of ore resources in 
mining and oil deposits in oil exploration, pollution evaluation in environmental 
sciences, to hydrology, meteorology, agriculture, etc. [H [2j [3] . For examples, en- 
vironmental resource models carry out spatial statistical analysis in the quantity 
of resources available such as the volume of available water, forest, etc., or their 
quality such as concentration of contaminants in air, water or soil samples. Ran- 
dom fields and their covariance functions or equivalently their variograms are used 
widely in the modeling of observed spatial data as these data are likely to be spa- 
tially dependent. The earlier developments of the subject include work by Whittle 
[HH], Matern Tatarski [8], Matheron [9] and others. The Gaussian random 

fields defined using the covariance functions from the Whittle-Matern (WM) co- 
variance class are widely used to model isotropic spatial processes in two and three 
dimensions. 

The WM class of covariance functions [TUl [TIJ H2 E2 03] has recently received 
considerable interest in geostatistics due to its great flexibility for modeling the 
spatial variations, in particular its ability to model behaviors of empirical variogram 
near the origin. Unlike other popular covariance models, the WM model has a 
parameter that characterizes the smoothness of the associated random field. Due 
to this reason, Stein [15] strongly recommended the WM class for the modeling of 
spatial covariance. 

A special case of the WM model was first obtained by Whittle [4] , who showed 
that a Gaussian random field with covariance function belonging to WM class can 
be obtained as a solution to a stochastic differential equation. The general form of 
WM model was given by Matern [6] and Tatarski [8], and was also considered by 
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Matheron [9] and Shkarofsky 16J. It can be associated with von Karman spectrum 
[T71 [18] in the modeling of wind speed. A comprehensive historical account on the 
WM class was given by Guttorp and Gneiting [T5] , who first called such a class of 
covariance functions as WM covariance family, but they later changed it to Matern 
covariance family [20 . 

Recall that the smoothness of a random field is characterized by the fractal di- 
mension, a local property which is determined by the asymptotic properties of the 
covariance near zero lag and its value depends on the smoothness parameter of the 
WM covariance. On the other hand, the strength of the spatial correlation is de- 
termined by the scale parameter and for large time lag it decays exponentially In 
this paper, we proposed a new generalization of WM covariance class with an addi- 
tional parameter which plays the role of scale or memory parameter, and the spatial 
correlation strength for large time lag now varies hyperbolically, with exponential 
decay as a special case. 

In the next section we recall some basic facts on the WM covariance class and 
the random field associated with it. This will be followed by the introduction of 
the GWM covariance class and the corresponding random field. The asymptotic 
properties of the GWM covariance function are studied in section [3l Based on 
these properties we are able to obtain the fractal dimension of the graph of the 
random field in GWM model. This random field satisfies a weaker self-similar 
property called local self-similarity, and it is short-range dependent. Simulations 
of the GWM covariance function and the random field (in two dimensions) are 
given. In the subsequent section, GWM process is applied to model wind speed 
and compared to model provided by WM process. Other possible generalizations 
and applications of the GWM random field are discussed in the concluding section. 

2. Generalized Whittle-Matern Model 
The WM class of covariance functions is given by [U El El |T5] : 

ID c m . (U) K^m, 

where K v {z) is the modified Bessel function of second kind (or Macdonald function), 
t £ K™, \t\ = + . . . + t\ is the Euclidean norm of t, A > is a scale parameter 
controlling the spatial range of the covariance, and v = 7 — (n/2) > is the 
smoothness parameter governing the level of smoothness of the associated Gaussian 
random field Y(t). Note that the WM covariance ((T|) has the same functional form 
as the characteristic function of the multivariate ^-distribution 2 lj . The spectral 
density of Y(t) is given by the Fourier transform of (fT|): 

The Gaussian random field Y(t) with covariance |T]) can be obtained as the solution 
to the following fractional stochastic differential equation 0]: 

(2) (—A + A 2 ) ^ Y(t) = rj(t), 

where A = Jp- + . . . + Jp- is the n-dimensional Laplacian, and r](t) is the standard 
white noise defined by 

(3) (*?(*)> = 0, (r)(t)r)(s)) = 5(t-s). 
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One early generalization of WM family of covariance functions was proposed by 
Shkarofsky |16j . Based on the argument that a covariance function for turbulence 
needs to have no cusp, it is required to have zero derivative at the origin and a 
second derivative that is finite and negative. In order to satisfy these requirements, 
he generalized ([l} to a covariance with two complementary parameters: 



(4) c(t) = «)<<*-,«) 

Clearly, up to constants, (HJ reduces to the covariance in WM class ((T|) when £ — > + . 
There also exist generalizations of WM class to a non-stationary class of covariance 
functions that allows for anisotropy, one such generalization is [22j : 



(5) C(ti,t 2 ) = 



A(ti+*a) 




In view of the wide applications of fractal operators in physics [23] . we pro- 
pose another generalization of WM class of covariance function by extending the 
fractional stochastic differential equation ([2]) to one with two fractional orders: 

(6) [(-A) a + Af r a , 7 (<) = #, 

with A, 7 > and a E (0, 1], and the Riesz fractional derivative D 2a = (— A) Q is 
defined by: 

(7) D 2a / = (-A)"/ = F- 1 {\u\ 2a F[f](u)} 
or 

(8) (FD 2 «/)(u,) = \u;\ 2a (F{f})(u), 

One can regard the fractional operator [(—A)" + A 2 ] 2 as a "shifted" Riesz deriv- 
ative and it has formally the series representation: 

(9) [(-A) a + A 2 ]^f; ( ll j 2 ) A) * 



See reference [26] for a more rigorous treatment of this operator based on hyper- 
singular integrals. Now by using 

(10) ([(-A) Q + A 2 ] * /) (t) = F- 1 ([|uf* + A 2 ] * F[/](uO) (t), 
the solution to ^ is found to be 

(11) Y ai {t) = ——^ / ^^d n u;, 

V ' M ' (2*)* J (|u;| 2 « + A 2 )* 

]SL n 

where f](u)) = F[i]](ijj) is the Fourier transform of the white noise. For convenience, 
we call Y an (t) the GWM (generalized Whittle-Matern) field. The representation 
(fTTj) shows that the GWM field Y an (t) is a centered Gaussian field with covariance 
function 

(12) C a , 1 (t-t') = C a , 7 (t,t') = (Y a , 1 (t)Y a ,,(t')) = ^ J dJp+^f ^' 
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From this, we see that Y a ^(t) is an isotropic field with spectral density 
(13) Sa ^ ] = (27r)«(M 2Q + A 2 ) 7 - 

Note that in the case of a = 1, 7 > 0, the field Y\ n (t) is sometimes called 
Bessel field by some authors [27J |28j [29j [30] based on the fact that the opera- 
tor (—A + A 2 ) 7 ^" is closely related to the Bessel potential with the WM covariance 
([l} equals to the Bessel kernel up to a multiplication constant [24] . Since the func- 
tion (|u;| 2q + X 2 Y l/2 is in L 2 (M) if and only if a 7 > n/2, the field Y an (t) is only 
well-defined by as an ordinary random field when > n/2. When a"/ < n/2, 
Y a ,-y(t) can be regarded as a generalized random field over the Schwarz space of 
test functions |31j . In the following, when we study the properties of the random 
field Y an (t), we restrict to the case > n/2. The two-dimensional GWM field 
with selected values of a and 7 are simulated in Figure [TJ In next section we shall 
study the asymptotic properties of the covariance and the sample path properties 
"fV„. .(/;. 

Here we would like to remark that when n = 1 and a = 1, the GWM process 
Yi )7 (t) is also called the Weyl fractional Ornstein-Uhlenbeck process or the Weyl 
fractional oscillator process [32], [33] , which can be considered as generalization of 
ordinary oscillator process driven by white noise. 
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3. Asymptotic Properties of the Covariance Function 

When a = 1, the covariance of GWM field (JT2J) Ci, 7 (t) reduces to the WM class 
given by |T]). However, (Til?]) in general does not have a closed analytic form. It is 
interesting to note that the spectral density of the GWM field has the same func- 
tional form as both the characteristic function of generalized multivariate Linnik 
distribution [34j [35] and the covariance function of generalized Cauchy class in M. d 
[36l[37]. Thus the covariance of the GWM field, the generalized multivariate Linnik 
distribution, and the spectral density of the random field belonging to the gener- 
alized Cauchy class all should have the same analytic and asymptotic properties. 
These properties have been considered for the generalized Linnik distribution in R, 
and the multivariate Linnik distribution for the special case with a G (0, 1) and 
7=1, and for the spectral density of the random field of generalized Cauchy class 
by Kotz et al. [33] and Ostrovskii [35], and Lim and Teo [37J respectively. Thus 
the results obtained in [331 (351 (37J can be translated directly to the covariance of 
the GWM field. 

For general a and 7, we can use a theorem of Bochner 38: which says that being 
an isotropic covariance function, C a , 7 (t) has a spectral representation given by 



(14) 



f°0 Jn-2 (Uj\t\) 

C a „(t) = (2tt)* / t ' , _ a Sa^U^dLJ 

Jo Mt|) 2 

f°° J^iHM) „ 

> 2 duj. 



(27r)fi (w 2q + A 2 )t 



Here J v (z) is the Bessel function of the first kind of order v. Now the result on 
a representation of the spectral density of the random field of generalized Cauchy 
class in [37] can be applied and we find that for a E (0, 1), the covariance function 
C aj7 (i) has another representation given by 



[tig f°° Kn^(u\t\) 

Jo W a « 2Q + A 2 ) 7 ' 



(15) C a , y (t) = - ^ Im / ——^——^u.du. 



In fact, for all a € (0, 1) and 7 > 0, the integral in (fTS"]) is convergent when t 7^ 0. 
Together with (fTJ), we find that for 0:7 < n/2, Ya i7 (t) can be considered as a random 
field with infinite variance and with covariance given by (|15[) if a G (0, 1); and by 
{TJ) if a = 1. Since 

/ 7T _ 

K v (z) ~ 4/ — e 2 as z — > 00, 

V 2z 

we can use (fTS"]) to effectively calculate the numerical values of C Qi7 (t). On the 
other hand, we can also use ([15]) to study the large \t\ behavior of the covariance 
function C an (t) when a € (0, 1). More precisely, using the formula 



E 



r(7 + j)(-i)-' j 



( i + z) 7 r( 7 ) i! 



and the formula 



f ,^,, )rf ^ 2 -r(i±f^)r(i±£^ 
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Figure 2: The graph of AC an (\t\), A = 2 2a, 7 r(a + l)T(a + n/2) sin(7ra)/(7r (n+2)/2 A 27+1 ), 
as a function of \t\. The reference curve is y = Itl" 2 " - ™. 



([39]. #6.561, no. 16), we find that if a € (0, 1), then when \t\ — > oo, we have 
(16) 

\t\~ U f°° K-n-2 (u) 



r^2 n+2 
2 2 7T 2 



Im 



o (e-yf^ + A 2 ) 



FiW'-tg 1X7) 



r(7 + j) (_J)l A -2 7 -j e i ff aj| t |-2aj /" Jjf^ ( u )d- 



V E r( S JH fa + !) r + I) 2^A-^ sin(^)|*r 2 ^-". 



3 t! r (^) J 
In particular, the leading term of C a ~(t) when \t\ — > oo is 

02a\-2f-l . 

(17) C a , 7 (t) J -T(a + l)r (a + - ) sin(7ra)|t|- 2Q - n . 

7T 2 \ Z / 

Note that the order of the leading term |t| -2Q ~™ only depends on a. In other words, 
the large time asymptotic behavior of the covariance function varies as |i|~ 2a ~™ and 
does not depend on 7. When a = 1, we cannot use (I15[) . However, we can obtain 
the large- 1 1 1 behavior of Ci )7 (i) from the explicit formula (Q} and the asymptotic 
formula for K v (z) ([32], #8.451, no. 6) which give 

(18) 



7T 2 1 (7) 



|t|-J+7-^ 



Notice that in this case, C\ 7 (i) decays exponentially and the leading term is 



(19) 



Ci, 7 (t) 



tt 2 ^ - r( 7 ) 



--^1*1 1 + 17- 



To study the local properties of the GWM field Y an (t) such as Holder continuity, 
local asymptotic self similarity and Hausdorff dimension of the graph, we need to 
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study the small- |t| behavior of the variogram 

(20) < 7 (i) := ([y„ l7 (t) - ^, 7 (0)] 2 > 
of the increment field Y a j (t) — Y a „(0). Notice that 

(21) < 7 (*) = 2(C a , T (0) - C a , 7 (t)), 
and the variance of Y^, )7 (t) is given explicitly by 

'[r«, 7 (t)] a ) =c a , 7 (o) = — L- / _J__ t r , « 



(22) = i r 

A g- 27 r ( 7 -£)r(£) 

2™7r*ar(f) r( 7 ) 

The |£| — > asymptotic properties of ct 2 7 (t) depend on the arithmetic nature of a 
and 7. To explore the leading behavior of cr 2 7 (t) as \t\ — > 0, we have to discuss the 
cases 07 e (§, I 4 :2 -), 07 = 2i j 2 - and cry > rj 4^- separately. From (J3TJ) and (TH| , 

2 -2 Z 100 / Jn-2 (fc|t|) 1 \ 

(23) ^{t) = ^j Q [ {k 2 lt{) ^ - 2 ^ r(f) J (fc2a+A 2 )7 *' 

Case I. When 0:7 € (§, ^y^), ^y making a change of variable fc 1— > (|2"5| is 

transformed to 



(24) (27r)f io ^ fc V 2^r(f)J (k^ + \z\t\^r dk - 

When \t\ — > 0, the integral 

/•«= / J„-2 (k) l \ fcn-l 

(25) m = Jo \ J^TWFV dk 

approaches a finite limit given by 

Using regularization method (see appendix), it can be shown that 

r(f-a 7 ) 

^ ' ~~ 2 2«7-fr(a7)' 

Therefore, as \t\ — > 0, 

( 27 ) = - 22a7 -i^ r( ^7 ) |t|aa7 " n + °(i*i 2Q7 "")- 

Notice that when a € (§ , ^y 2 -), the leading order of a 2 (£) depends on a and 7 
only in the combination cry. By letting 7 = 7 '/a gives 07 = 7'. Hence the \t\ — > 
asymptotic properties of cr 2 (t) vary as \t\ 21 ~ n which is independent of a. 
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Figure 3: The graph of Bcr 2 a ^(t), B = -T(n/2 - aj)/(2 2a ' 1 ' 1 TT n/2 r{aj)), as a function 
of \t\. The reference curve is y = \t\ 2aj ~ n . 



Case II. When > ^t^, using the fact that 

z 



J ' (z) = ^^^7!r>T7Ti)- 



([33, #8.402), we find that as |t| 0, 
(28 ) 1 



o{\t\ 2 ). 



(29) 



2^r(f) 2^r(^) 

Therefore, (|2"3")1 gives 

2 , , _ |t| 2 ^ fc" +1 rffc 2 

^W-^fr^) X (fc 2a + A 2)7 ) 

"2«+i7Ttar(^) r( 7 ) 11 + u 1 j 

as |t| -> 0. 

Case III. In the limiting case 07 = ^y 2 -, eq. (|24|) gives 



(30) o* (t) 

However, now the integral 

r°° / J„_ 2 (fc) 
J(t) = ' — 



-2|t| 2 Z" 00 / ^(A) 1 



1.71— 1 



-dk. 



1 



n — 1 



-rife 



/o V k ^ 2 ^ r (f) / (fc 2Q + A 2 |«| 2a r 
does not have a finite limit when t — > 0. In the appendix, we show that 



(31) 



for some constant A. Therefore, as \t\ — > 0, 
(32) - 2 



1 .-2, 1 



^ (t) = 2 >fr(f) |f|2log W +Q(|f| 
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From |27|) . (|29|) and (|32|) . we see that the behavior of the leading order term of 
cr^ 7 (t) when |t| -> depends on 7' - (n/2). If 7' - (n/2) <E (0,1), the leading 
order term is of order |i| 27 ~ n which depends on the magnitude of 7' — {n/2). If 
7' — (n/2) > 1, then the leading order term is of order \t\ 2 , which loses dependence 
on 7' — (n/2). In the borderline case 7' — (n/2) = 1, the leading order term is of 
order |*| 2 log(l/|*|). 

The graphs of C a , 7 (|i|) when \t\ is large and cr 2 7 (|i|) when \t\ is small for some 
particular values of a and 7 are given in Figures [2] and [3] respectively. 

4. Sample Path Properties of GWM Field Y an (t) 

Some basic sample path properties of the GWM field will be considered in this 
section. 

4.1. Continuity and differentiability. When > n/2 for which the field 
Y a ^(t) is defined as an ordinary random field, the covariance function C Q;7 (t) 
(TT2"]) is continuous at t = 0. By a well-known result (see e.g. [40J ) . this implies 
that the field Y a ^(t) is mean square (m.s.) continuous. One may then proceed 
to investigate the differentiability of the field Y an (t). It turns out that Y an (t) is 
not always differentiable. In fact, a well-known result (see e.g. [40]) states that the 
m.s. first partial derivative dX(t)/dtj of a stationary random field X(t) exists if 
and only if the partial derivative d 2 C(t)/dt 2 exists at t = 0, where C(t) denotes 
the covariance function of X(t). From our result in the previous section, we find 
that as t — > 0, 

C(t) - C(0) 



(33) 



'B 1 |t| 2 o7-« + (|t| 2a T- n ), if a 7 e (f 



. 2 ' 2 



= <J J B 2 |<| 2 log^+0(|t| 2 ), ifa 7 =^, 
B 3 |t| 2 + o(|t| 2 ), ifa 7 >^, 



for some constants B\, B2, B3. It is easy to check that for the radial function 
f(t) = \t\ , the second partial derivative d 2 f(t)/dt 2 exists at t = if and only if 
h > 2. Therefore we conclude that the mean square partial derivatives of Y an (t) 
exist if and only if aj > , with a representation given by 

dY an i f ujje^^uj) 



(34) ^^(t) = -— / ^^(Ta;. 

dtj (2ir)~ J (| w |2a + A 2)2 

In fact, we can argue analogously that the m.s. j-th order partial derivatives of 
Y an (t) exist if and only if > (n/2) + j. 

For «7 € (s, ^i 2 -) , the field y Qj7 (t) is not differentiable. Therefore, we would 
instead investigate the order of continuity of Y" a , 7 (t). Recall that a function / is 
said to be Holder continuous of order h G (0, 1] if and only if 

(35) \f(t')-f(t)\<K\t'-t\ h V t',t 

for some constant K. The sup of all h where / is Holder continuous of order h is 
called the Holder exponent of /. For a centered isotropic Gaussian random field 
X(t), a concept of index-/? field was introduced by Adler [40] which can be used 
to characterize the Holder exponent of the sample paths of X(t). More precisely, a 
theorem states that if X(t) is an index-/? field, then with probability one, its sample 
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paths have Holder exponent equal to j3, where X(t) is called index-/? field if and 
only if 

/?=sup(/3 : cr(t) = o(\tf) as \t\ -> o) 

(36) r- a i 

= inf|/3 : Itl' 3 = o(<t(*)) as |t| -► 0| . 

Here er(£) is defined as the square root of the variogram of X(t) 1 i.e., a(t) = 

[X(t) - X(0)] 2 y For the field F a , 7 (t) we are considering, it is immediate to 

conclude that from (H7J), @ and ((32J) that if 07 G (f, 1 ^), then Y an (t) is 
an indexed (017 — [n/2)) field; whereas if 07 > then Y" Q . 7 (i) is an index-1 

field. Therefore, we have for cry £ (^, the sample paths of Y an (t) is Holder 

continuous of order 07— (n/2) with probability one. For > ^y^, it can be shown 
by considering the gradient field \7Y an (t) — (dY a ^(t) / dt\, . . . , dY an (t) /dt n ) that 
the sample paths of Y an (t) are diffcrcntiable. 

4.2. Fractal dimension. For a non-differentiable function /, ordinary definition 
of dimension, which is always a nonncgative integer, is inadequate to measure the 
dimensionality of its image or graph. A more appropriate definition of dimension is 
called fractal of Hausdorff dimension which can be any non- negative real number. 
Definition and basic properties of fractal dimension can be obtained in the book 
[41 j . Here we would like to make use of the following result. For an index-/? field 
in K™, with probability one, the fractal dimension of the image and graph of its 
sample path are 1 and n+ 1 — (3 respectively. Thus, with probability one, the image 
of the sample path of Y a ,-y(t) always has fractal dimension one. The result is more 
interesting for the fractal dimension of the graphs. If aj G (§, ^2 )' then with 
probability one, the graph of the sample path of Y a>7 (t) has dimension + 1 — cry, 
a real number between n and n + 1. However, when 07 exceeds the point ^4^, 
then with probability one, the graph of the sample path of Y an (t) always has 
dimension equal to n. This is reasonable since the sample path of Y Q)7 (i) becomes 
differentiable when aj > ^J^- In fact, Figure [1] show clearly that the fractal 
dimension of the graph of Y an (ti, t%) depend on 7' = cry. 

4.3. Local self-similarity. Self-similarity is an important property of fractals. 
Intuitively, a field is called self-similar if it is invariant under appropriate scaling. 
For a random field X(t), we say that it is self-similar of order H if and only if for 
any r > 0, the law of the field X(rt) is the same as the law of the field r H X(t). It is 
well-known that a stationary random field cannot be self-similar [42] . In fact, up to 
a constant multiplicative factor, the only ii-self-similar centered Gaussian random 
field with stationary increments is the fractional Levy Brownian field Bh (i) of index 
H with covariance 

(37) (B H {s)B H {t)) = \ {\t\ 2H + \s\ 2H -\t- s\ 2H ) . 

This excludes the possibility for Y an (t) being a self-similar random field. How- 
ever, Y an (t) satisfies a weaker self-similar property known as local self-similarity 
considered by Kent and Wood A centered stationary Gaussian field is locally 
self-similar of order (3/2 if its covariance C(t) satisfies for |t| — > 0, 

(38) C(t) = C(0)-A\tf[l + O(\t\ s )] 



GENERALIZED WHITTLE-MATERN RANDOM FIELD 



11 



with A > and 5 > 0. The proof of eq. (27]) shows that for E (§, 
C a , 7 (t) = C7 O)7 (0) - A\t\ 2a ^ n + o(|t| 2 ^-"+ 5 ) 

with 

(39 ) A = 1 F (t-"7) 

Hence Y an (t) is locally self-similar of order aj — (n/2). 

There exists an equivalent way of characterizing self-similarity at a local scale 
called local asymptotical self-similarity which was first introduced for multifrac- 
tional Brownian motion [33]. Recall that a random field X(t) is called locally 
asymptotically self-similar with parameter H G (0, 1) at a point to if the limit 
random field 

f X(t + pu) - X(t ) 

(40) |T i0 («) = \im + K0 , u e 

exists and is nontrivial [33]- In this case, Tt (u) is called the tangent field of X(t) 
at to. It can be directly verified that 

([Y an (to + pu) - Y an (t )} [Y a>1 (t + pv) - F a , 7 (t )]) 

(41) ; 

2 

Eq. (H7J) then shows that for aj € (- ^ 



= o « 7 (/> M ) + °^, 7 (H " < 7 M M ~ «))) 



lim 



2 ' 2 / ' 

F(i + pu) - Y(t Q ) Y(to + pv) - Y{t ) 



(42) p-o+ \ p a ~>~2 

=A (\u\ 2a ~<- n + |«| 2Q T- n - | u - v\ 2a ~<- n ) . 

Therefore, for ct"f € (fi 2 ^" 2 ")' Y an (t) is locally asymptotically self-similar with 
order aj — (n/2). Its tangent field at a point to € M" is independent of to, and up 
to the multiplicative factor 2A, the tangent is given by the fractional Levy Brownian 
field B H (u) d3Tj) of order - (n/2). 

From the results on Holder exponent, fractal dimension and local self-similarity, 
it is found that they all depend on the parameters a and 7 in the combination 
7' = cry. They are related to each other in such a way that if the Holder exponent 
is H = min{7' — (n/2), 1}, then the fractal dimension is n + 1 — H and the order 
of local self-similarity is H again if H < 1. 

4.4. Short Range Dependence. Recall that a stationary random field X(t) is 
said to have short range dependence (or short memory) if the absolute value of its 
covariance function C(t) is integrable over K™, that is 



J \C(t)\d n t < 00. 



By our result on the large \t\ asymptotic behavior of C a ^(t) (fT7|) and (fT9| . we find 
that when \t\ -> 00, C an (t) ~ \t\~ 2a ~ n if a G (0,1) and C an (t) ~ e-^t^-^ 
if a = 1. Using polar coordinates, it can be verified easily that 

\t\ p d n t < 00 

tem™,|t|>i 
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if and only if p < —n. This immediately implies that Y a ^(t) has short range de- 
pendence for all a and 7. Moreover, for a £ (0,1), the short memory exponent 
2a + n depends only on a and not on 7. Together with the result on local properties 
such as Holder exponent, fractal dimension and local asymptotic self-similarity, this 
implies that the short range dependence property and local properties of the GWM 
field Y a ,j(t) are characterized separately by a and 7' = cry. This should be com- 
pared to the random field with generalized Cauchy covariance [36, 37J, which is an 
isotropic random field with two parameters that enables separate characterizations 
of long range dependence and fractal dimension. 

Here we would like to remark on the Markov property for WM field and GWM 
field. In the case of Whittle field (a = 7 = 1) in n = 2 dimension, this problem 
has been studied by Pitt and Robeva (45] , who showed that under certain technical 
conditions, the sharp Markov property is satisfied. They generalized the result 
to WM field (which they called Bessel field), and they verified that under some 
technical conditions the sharp Markov property holds for WM field with n + 1/2 < 
7 < n + 1, n > 1 28, 29]. It will be interesting to see whether the arguments of 
Pitt and Robeva can be extended to the GWM field. 



WM field has been widely used in modeling [TUl HH 021 Q21 [14] geostastical data 
such as sea beam data, temperature, wind speed and soil data. In this section, we 
show that the GWM process can be used to provide an alternative model for wind 
speed. 

We analyze the average daily wind speed of Roche's Point in Ireland from 1973 
to 1978 which consists of N = 365 x 6 = 2190 data pointfQ (see Figure HJ). The 
Irish wind data of 12 meteorological sites from 1961 to 1978 has been analyzed by 
several authors [321 H7] HH1 H2] where they were more concerned with the spatial 
correlation between the sites. On the other hand, the von Karman wind turbulence 
model [T71 [18] proposed a model with spectral density having the same functional 
form as the spectral density of the WM process. 

As in [46] , we consider the seasonal effect by calculating the average of the square 
roots of the daily means over the 6 years for each day of the year, and then regressing 
the result with a polynomial of degree 8 (see Figure . The deseasonalized data 
(Figure [6]) is obtained by subtracting the fitted polynomial from the square roots 
of daily means. It has zero mean and is referred to as the velocity measures. To 
justify that the velocity measures is a short memory process, we use the fact that 
a process has long memory if and only if it's spectral density diverges at u = 0. 
For a discrete stationary random process X tl t — 1, 2, 3, . . ., with covariance C(t), 
an analog of spectral density is the power spectral density (PSD) defined by 



It is a periodic function with period 2-7r and S(2n — u>) — S(lo). If X t has an 
underlying continuous process X(t) with spectral density S(uj) so that X t — X(t) 



5. Application to wind speed modeling 




^The data is obtained from Statlib (http:/ /lib. stat.cmu.edu/datasets/ ) with the value for 29th, 
February, 1976 omitted. 
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Daily average wind speed at Roche's Point, Ireland 

20 1 1 — i 1 1 1 




QJ , , , , , 

1973 1974 1975 1976 1977 1978 1979 



year 

Figure 4: The daily average wind speed at Roche's Point, Ireland from 1973 to 1978. 




Figure 5: The average of the square roots of daily means over the 6 years for each day of 
the year and the fitted polynomial of degree 8. 



when t — 1,2,3,..., then 



PSD(w) = J2 S ( UJ - 27r -?)- 



J = -CO 



There are different ways to estimate the power spectral density from a given sample 
x t ,t = 1, 2, . . . , N of X t . One way is to use the periodogram method, where the 
estimate of PSD(w) is given by the periodogram 



1 

27riV 



JV 
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year 



Figure 6: The deseasonalized data (velocity measures). 




Figure 7: The estimated power spectral density of the velocity measures and the corre- 
sponding log-log-plot using (a) periodogram method, (b) Welch's method. 



However, this method usually leads to large fluctuations. A better method which 
gives a smoother estimate is introduced by Welch [50 and improved by others (see 
e.g. [51]). We estimate the PSD of the velocity measures using Welch's method 
by segmenting the data into 50% overlapping blocks of length 73 and applying the 
Hamming window to each block. The resulting estimate for PSD is compared to 
the periodogram estimate in Figure [7] The behavior of PSD at u> ~ gives strong 
evidence that the velocity measures has short memory. 

In order to apply the GWM process to model the velocity measures, we consider 
the more general process parametrized by four parameters a, 7, K, I: 



Y£f(t) = KY a ^(£t) x=1 = -= / - " ' da;, 

/ ■) — t .r'.rv 1 ^ \ 2 



which has spectral density 
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and covariancc function 



C^(t) = K 2 C an (£t) x= 



Notice that I rescale the time parameter and K rescale the data. We need to 
determine the parameters a,^,K,£ so that the process Y^ e (t) gives the best 
model to the velocity measures. For this purpose we use the maximum likeli- 
hood estimation (MLE) method strongly recommended by Stein (see e.g. [15]). Let 
T(0) = r(a,7, K, £) be the covariance matrix (C^(i — j))^ =1 . Since we assume 
that the velocity measure is a Gaussian process, the probability density function 
for y = (yi, . . . , un) t having mean and covariance Y{ff) is 

(43) p(y,0) = 1 Pyp ( -lyTrmy^j 

V ; V ' (2tt)£ v/ditr(0) V 2 y v ' 

In MLE method, we seek the parameters 9 — (a,j,K,£) that would maximize 
the probability density function (|4"3"| with y being the observed velocity measures. 
Equivalently, we have to minimize the negative log of the likelihood function: 

(44) NLL(O) = - logp(y; 9) = ^y T T(0y 1 y + ~ logdct T(0) + y log(27r). 

Finding the minimum of the highly nonlinear function (144)) with four parameters 
is computationally demanding. Therefore it is desirable to reduce the number of 
parameters which has to be estimated. Notice that the variance s 2 of Y£~ (t) is 
given by 

,2_nKi^_ k 2 r(i)r( 7 -i) 



* 2 = «(0) 



2na 

As a result, the value of K can be determined from this equation once s 2 ,a,7 
are given. On the other hand, we can rewrite the covariance matrix T(6) — 
T(a,-f,K,£) as s 2 p(0') = s 2 p(a,j,£), where p(a,"f,£) is the correlation matrix 

(p%(i-i))Lv 



which is independent of K. Rewriting in the variables a,j,£, s 2 , we have 

(45) NLL(G', s 2 ) ^y T p(6 , y 1 y + ^ logs 2 + \ logdet p{6) + ^ log(27r). 



Taking derivative with respect to s 2 , we find that for fixed 6' = (a, "f,£), the 
minimum of NLL(9\ s 2 ) appears at 

(46) s 2 = ±y T P {e')- l y. 

Substituting this into (|45|) . we reduce the problem to finding 0' to minimize the 

function 

(47) 

NLL{6') =j logy T p(ey 1 y + ± logdct p(9') + y (1 + log(27r) - log(iV)) , 

and s 2 is then determined from (|46[) . The fminsearch function in Matlab which 
uses the simplex search algorithm by Neldon and Mead [S3] , is used to identify the 
minimum of (|47[) . This algorithm does not involve computation of derivatives. In 
order to compare the GWM model with that of WM, we also run the same search 
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Figure 8: The empirical power spectral density of the velocity measures compared to the 
PSD of the WM model and the GWM model with parameters given in Table [TJ and the 
corresponding log- log-plot. 



Table 1: The estimated parameters for the WM and GWM models 





a 


7 


K 


t 


s 2 


NLL 


WM model 


1 


1.0225 


0.7857 


0.7474 


0.2994 


1488.42 


GWM model 


0.5186 


4.1223 


1.6965 


2.8250 


0.2995 


1487.47 



with a fix to 1. The results are tabulated in Table [TJ It suggests the WM model 
Y^' (t) with spectral density 

0.5796 1 



2tt (lu 2 + 0.75 2 ) 1 02 
and the GWM model Y^ e (t) with spectral density 

50.9376 1 

AGWMM - 2n (M 1.03 + 2.821.03)4.12 

for the velocity measures. From Table [TJ we see that the GWM model gives a 
better value to NLL(0). On the other hand, a graphical comparison of the PSD 
of the WM model and the GWM model for velocity measures and the empirical 
PSD (Figure [8j) also shows that the GWM model gives a better fit to the velocity 
measures compared to the WM model especially in the low frequency region. 
Here we would also like to remark that theoretically, the variogram cr^(/i) 2 = 

Y£f(t + h)- Y a K /(t)] 2 ) = 2{C«m - C^Mh)) approaches 



^ (0) _^ 2 r(^)r( 7 -^; 

Q ' 7 ira r(7) 
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0.45 
0.4 
0.35 



500 1000 1500 

h 

Figure 9: The empirical variogram a 2 (h). The horizontal line gives the estimated value 
of C*'y(0) = 0.2964. 




10° io' io 2 

h 



Figure 10: The empirical variogram 5" 2 (/i) compared to the variograms of the WM and 
GWM models. 



as h — > oo. Figure [9] shows the empirical variogram of the velocity measures esti- 
mated by 

N-h 

v 2 {h) = N _ h ( yi+h ~ Vi f ■ 



The horizontal line gives an estimation of the variance s 2 = 0.2964, which is very 
close to the one estimated by MLE. Figure ITOl compares the empirical variogram to 
the variograms of the WM model and GWM model for velocity measures. 
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6. Concluding Remarks 

In this paper, we have introduced a new class of Gaussian random field with co- 
variance belonging to a generalized Whittle-Matern family of covariance functions. 
Some of the basic properties of this GWM field are studied. Simulations of the 
GWM covariance and GWM field in two dimensions are carried out. We also apply 
this random process to model wind speed. It is shown that this new random field 
can provide a more flexible alternative to modeling. In the future, we would like 
to extend the application of GWM field to provide models for other geostatistical 
data such as sea beam data, geothermal field temperature and soil data which the 
WM field has been shown to provide a good model [TQl HH H21 EH [M] . 

Just like its predecessor WM model, GWM model will find its main applications 
in geostatistics. However, one expects it can have potential applications in mod- 
eling short range dependent process such as coding regions of DNA sequences and 
fluctuations of an electropore of nano size [SSJ HH EZ1 HSJ IS9] 150] . It will also be in- 
teresting to consider its applications in modeling fractional diffusion and fractional 
anomalous diffusion [611 E2 O EH [65j [66j [67] . This later aspect is being consid- 
ered elsewhere. For d = 2, GWM may serve as a model to two-dimensional images 
arising in biology, chemistry and physics, in addition to those from geological and 
environmental images. The advances in imaging techniques allow better analysis 
of morphology of various material surfaces. Various spatial statistical and mor- 
phological methods are available to analyze the patterns of the surface of complex 
materials, hence the characterization of their physical properties. Examples of data 
that can be modeled by the random field include the concentration of particular 
component in a liquid or solid sample, properties such as porosity, permeability, 
conductivity, absorptivity, emissivity, etc. of the material samples. Although ap- 
plications of spatial models to statistical physics are still quite limited, some recent 
efforts have been made in [68] ESI HQ] • One expects GWM model also has such 
potential applications. In particular, as a correlation model, GWM model is use- 
ful in the modeling of morphological structure of complex material with spatial 
correlation that is short-ranged, that is, its underlying physical process is weakly 
correlated or weakly coupled over finite spatial or temporal scales. Readers can 
consult references in |71] for recent advances and applications of spatial models in 
physics, in particular statistical physics and astrophysics. 

For modeling of data which may have correlation with time or space dependent 
memory parameter or smoothness parameter, it is necessary to consider GWM field 
Y an (t) with a and 7 replaced by a(t) > and j(t) > 0, with a{t)^f (t) > n/2. 
A GWM field with two variable fractional indices can be studied in a similar way 
like the fractional Riesz-Bessel field with variable order [30]. In fact, when n = 1 
and a = 1, such a generalization has been considered in [53] where it is called 
Weyl multifractional Ornstein-Uhlenbeck process. Another possible extension is 
the anisotropic counterpart of the GWM field whose covariance is a product of 
GWM processes. 
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Appendix A. Derivations of formulas (26) and (31) 

1. We want to prove eq. (|25|) when c*7 £ (§, Using regularization method, 

we have 

J„^(fc)fcf l f oo f,n-l dk 

1 = hm < / -7TT. ^ — dk 



<*->0+l./ (fc 2 + a 2 )^ 2^r(f)7o (fc 2 + a 2 )" 7 ]' 

Applying the formulas #6.565, no. 4 and #3.251, no. 11 of [39], we find that 

a^o+ 1 2 Q T- 1 r(a7) 7 jl ' 2f r(a7) J 
Now the formulas #8.485 and #8.445 of [39 give 

2sin(7rz/) I ^j!r(j + ^j!r(j + l + i/) 

when j/ S (0, 1); which allows us to conclude that 

T 1 r(f-a 7 ) 



2 2 Q7 -f sin ( a7 _ l)] r (cry) r (cry - f + l) 2 2a ^-t r(a 7 ) ' 

2. We want to prove eq. (f3Tj) when 07 = Using ([28)1 . we can write I(t) as the 
sum of I\ (t) and I2 (t) , where 

r 1 f Jz^jk) 1 fc 2 \ fc"- 1 

,f, - / 1 ^ 2 ^r(f) + 2^r(^) y ) (fc 2Q + A 2 |t| 2 «r* 



afc 



fc^ 2^r(f); (fc 2a + A 2 |*| 2 «) 

has a finite limit /i(0) as \t\ — > 0, and 

1 Z" 1 fc" +1 
J 2 (t) := /(*) - ii(t) = T2 / 77^ ,„ ,0 s dk. 

2y ' K> n> 2^r(^)7o (fc 2 « + A 2 |*| 2 «)T 

By making a change of variable fc 1— > fc 1 ^ 2 "', we have 

fcT- 1 * 



2^aT(^)Jo (fc + A 2 |i| 2 «r' 



From this, we find that ^(t) can be written as a sum of /3(f) and I&(t), where 

1 r 1 f fc''- 1 1 1 

sW- ~ 2 ^ aT ^J Q \(fc + A 2 |t| 2 ^ ~ fc + A 2 |*| 2 «J 
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has a finite limit /3(C)) when \t\ — ► 0, and 



1 , 1 + A 2 |t| 2t * 

■log- 



2 ^ a r(^) 6 a 2 |*| 



2 a 



lQ g^ + 7^+4 j-7-T?T logA 2 + o(l). 



2 ^r(^) 1*1 2 2 F«r(if 

Therefore, we have shown that 

J ( f ) = - 2 ^ r 1 m 1 ° S M +A + 0(1) ' 

where 

* = AW) + 7,(0) + ^ kg*. 
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